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Abstract 

IT) ■ 

^v^j . In this paper, we address the problem of safety verification of nonlinear hybrid systems. 

A hybrid symbolic-numeric method is presented to compute exact inequality invariants of 
i i hybrid systems efficiently. Some numerical invariants of a hybrid system can be obtained 

w: 

by solving a bilinear SOS programming via PENBMI solver or iterative method, then the 
^0 1 modified Newton refinement and rational vector recovery techniques are applied to obtain 

, exact polynomial invariants with rational coefficients, which exactly satisfy the conditions 

O ■ of invariants. Experiments on some benchmarks are given to illustrate the efficiency of our 

algorithm. 

q^; 1. Introduction 

£Nl 1 Complex physical systems are systems in which the techniques of sensing, control, communication 

and coordination are involved and interacted with each other. Since many of such systems are 
7-H 1 safety critical systems, such as the controllers widely used in airplanes, railway, and automo- 

biles, ensuring correct functioning of these systems is among the most important and challenging 
problems in computer science, mathematics and engineering. As a common mathematical model 
for complex physical systems, hybrid systems [6] are dynamical systems that are governed by 
interacting discrete and continuous dynamics. Continuous dynamics is specified by differential 
equations, and for discrete transitions, the hybrid system changes state instantaneously and pos- 
sibly discontinuously. Among the most important research issues on hybrid systems are those 
of safety, i.e., deciding whether a given property holds in all the reachable states, and its dual 
problem reachability, i.e., deciding if there exists a trajectory starting from the initial set that 
reaches a state satisfying the given property. 

Due to the infinite number of possible states in continuous state spaces, safety verification 
or reachability analysis of hybrid systems presents a challenge. Some well-established techniques 
have been proposed. In [2], level set methods and flow-pipe approximations were presented for 
computing approximate reachable sets of hybrid systems. By contrast, quantifier elimination was 
used in [11] to compute exact reachable sets for linear systems with certain eigenstructures and 
semi-algebraic initial sets, and this method was generalized in [28] to handle linear systems with 
almost arbitrary eigenstructures. Recently, invariant generation has been proposed for safety 
verification of hybrid systems. An invariant [25] of a hybrid system is a property that holds 
in all the reachable states of the system, in other words, it is an over-approximation of all the 
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reachable states. Invariants are useful facts about the dynamics of a given system. However, 
generating invariants of arbitrary form is known to be computationally hard, intractable even for 
the simplest classes. The usual technique for generating invariants is to compute an inductive 
invariant, an assertion that holds initially and is preserved by all discrete and continuous state 
changes. There has been lots of work towards invariant generation for hybrid systems using convex 
optimization and semi-algebraic system solving [24, 5, 16, 19, 25, 20, 28, 27]. However, some of 
these techniques are subject to numerical errors and some suffer from high complexity. Taking 
advantage of the efficiency of numerical computation and the error-free property of symbolic 
computation, we proposed in [29] a hybrid symbolic-numeric method via sum of squares (SOS) 
relaxation and exact certificate to construct differential invariants for continuous dynamic systems, 
and generalized in [15] the idea for safety verification of hybrid systems. 

In this paper, we present a hybrid symbolic-numeric algorithm to compute exact invariants 
of hybrid systems. The algorithm is based on SOS relaxation [17] of a parametric polynomial 
optimization problem with bilinear matrix inequality (BMI) constraints, which can be solved 
directly using a recently developed PENBMI solver [9, 10] in Matlab or an iterative method, and 
exact SOS representation recovery techniques presented in [7, 8]. The algorithm in this work 
improves our former result in [15] in two aspects. First, we replace the strengthened linear matrix 
inequality (LMI) constraints with the original BMI constraints in the parametric optimization 
problem. Second, we modify both Newton iteration refinement and rational vector recovery, 
which can handle some cases where the method in [15] fails and usually yield invariants with 
lower degree. Unlike the numerical approaches, our method can yield exact invariants, which 
can overcome the unsoundness in verification of hybrid systems caused by numerical errors [18], 
as illustrated in Example 2. In comparison with some symbolic approaches based on qualifier 
elimination technique, our approach is more efficient and practical, because parametric polynomial 
optimization problem based on SOS relaxation can be solved in polynomial time theoretically. 

The rest of the paper is organized as follows. In Sections 2, we introduce some notions about 
hybrid systems and safety verification. And in Section 3, we transform the problem of safety ver- 
ification of hybrid systems into a parametric program with BMI constraints. Section 4 is devoted 
to computing numerical solutions of a BMI problem via PENBMI solver or iterative method. 
In Section 5, an algorithm based on the modified Newton iteration and rational vector recovery 
techniques is proposed to obtain exact solutions of BMI problem with rational coefficients. In 
Section 6, we discuss issues related to the implementation of the proposed mehtod. In Section 7, 
experiments on some benchmarks are shown to illustrate our algorithm on safety verification. 
Section 8 concludes the paper. 

2. Hybrid Systems and Safety Verification 

To model hybrid systems, we use the notion of hybrid automata [6, 25]. 

Definition 1 (Hybrid system). A hybrid system H : (V, L, T, O, T>, "J, £q) consists of the fol- 
lowing components: 

• V = {xi, x n }, a set of real-valued system variables. A state is an interpretation of V , 
assigning to each Xi G V a real value. An assertion is a first-order formula over V . A 
state s satisfies an assertion tp, written as s \= ip, if ip holds on the state s. We will also 
write tpi \= pi for two assertions <p\,pi to denote that <pi is true at least in all the states 
in which p\ is true; 

• L, a finite set of locations; 

• T, a set of (discrete) transitions. Each transition r : (£,£', g T , p T ) G T consists of a pre- 
location £ G L, a postlocation £' G L, the guard condition g T over V , and an assertion p T 
over V U V representing the next-state relation, where V = {x' x , ...,x' n } denotes the next- 
state variables. Note that the transition r can take place only if g T holds; 
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• ; an assertion specifying the initial condition; 

• T>, a map that associates each location £ G L to a differential rule (also known as a vector 
field ) 'D{£), an autonomous system xi = fi^{V) for each Xi G V , written briefly as x = f^(x). 
The differential rule at a location specifies how the system variables evolve in that location; 

• , a map that maps each location £ G L to a location condition (location invariant) ^(l), 
an assertion over V ; 

• £q G L, the initial location. We assume that the initial condition satisfies the location 
invariant at the initial location, that is, (= ^(^o)- 

Given a hybrid system H with a prespecified unsafe region X u C R™, we say that the system 
H is safe if all trajectories of H starting from any state in the initial set, can not evolve to X u , 
or, equivalently, any state in X u is not reachable. We can also specify an unsafe region, denoted 
as X u (£), for each location £ G L. 

For safety verification of hybrid systems, the notion of invariants of hybrid systems plays an 
important role. We define 

Definition 2 (Invariant). [25] An invariant of a hybrid system at location I is an assertion I 
such that for any reachable state (£, x) of the hybrid system, x (= X. 

An invariant of a hybrid system is an assertion that holds in all the reachable states of the 
system. 

Clearly, an invariant of a hybrid system is an over-approximation of all the reachable states of 
the system. If an invariant lies inside the safe regions, or its intersection with the unsafe regions 
is empty, then safety of hybrid systems is verified. However, generating invariants with arbitrary 
form is known to be computationally hard, intractable even for the simplest classes. The usual 
technique for generating invariants is to compute inductive invariants, as defined below. 

Definition 3 (Inductive invariant). An inductive assertion map I of a hybrid system H : 
(V, L, 7~, O, T>, ^, £q) is a map that associates with each location £ G L an assertion X(£) that holds 
initially and is preserved by all discrete transitions and continuous flows of H. More formally, 
an inductive assertion map satisfies the following requirements: 

[Initial] 6 \=1(£ )- 

[Discrete Consecution] For each discrete transition r : {£,£', g T ,p T ) starting from a state sat- 
isfying T{£), taking r leads to a state satisfying !(£'), i.e., I(£) Ag T A p T \= !(£') where !{£') 
represents the assertion T(£) with the current state variables X\, . . . , x n replaced by the next 
state variables x[ , . . . , x' n , respectively. 

[Continuous Consecution] For every location £ G L and states (£, xi), (£, X2) such that X2 
evolves from xi according to the differential rule T>{£) at £, if Xi |= T(£) then X2 |= !(£)■ 

By a polynomial hybrid system, we mean a hybrid system H : (V, L, 7, Q, T>, ^ , £0), where 
the initial condition 0, location invariants ^{£), and the guard condition and reset relation in 
each transition r G T are conjunctions of polynomial inequalities over the program variables, and 
moreover, each differential rule T>{£) is of the form Xi = /e,t( x ) with fi,i(x) G R[x]. 

In a preceding paper [15], we propose a symbolic-numeric approach to generate polynomial 
invariants of the form y(x) > for polynomial hybrid systems via the combination of Sum-of- 
Squares (SOS) relaxation with Gauss-Newton refinement and rational vector recovery techniques. 
We will describe how to improve our result in [15] by solving the BMI problem directly. As stated 
in the following theorem, safety verification of hybrid systems can be reduced to finding invariants 
of hybrid systems. 

Theorem 1. [Theorem 2.7, [20]] Let H : (V, L, T, Q,T), ^,£0) be a hybrid system. Suppose for 
each location £ G L, there exists a function ^(x) such that the following conditions hold: 
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(i) e h^o(x) >o, 

(ii) </?«( x ) > A g(£,£') A p(£,£') (= </^'(x') > 0, /or any transition (£,£', g, p) going out of £, 

(iii) <^(x) = A $(1) |= ^(x) > 0, here ^(x) denotes the Lie-derivative of <Pi(x) along the 
vector field T>(£), i.e., <pt(x) = Ya = i §ff ' /«,i( x ); 

iften V^( x ) > is an invariant of the hybrid system H at location £. If moreover, 

(iv) X„«) h WW < 0, We L, 

i/ien i/ie safety of the system H is guaranteed. 

Proof. The proof is obvious. □ 

Clearly, to guarantee safety of the hybrid system H, the intersection of the initial set O and 
the unsafe region X u (£ ) must be empty. When the hybrid system H in Theorem 1 specializes 
to a continuous system at one location, denoted as D : (V, 6,1?,*), then the condition (iv) in 
Theorem 1 can be relaxed through replacing the whole unsafe region X u by its boundary dX u , 
as illustrated by the following. 

Corollary 1. Let D : (V, 0,2?,*) be a continuous system. Suppose there exists a function <£>(x) 
satisfying the following conditions: 

(i) e h v(x) > o, 

(ii) 9?(x) = A * h <p(x) > 0, 

(iii) dX u |= </?(x) < 0, Ziere dX u denotes the boundary of the set X u , 
then the safety of the system D is guaranteed. 

Proof. By conditions (i) and (ii), the values of y(x) can not be negative during the entire evolution 
of the system D. Then condition (iii) implies that all reachable sets lie outside the unsafe region 
X u , yielding the safety of the system. □ 

In addition, we can consider more complicated forms of invariants for safety verification, i.e., 
invariants that are conjunctions of several polynomial inequalities: /\ t (pe t i(x) > 0. For simplicity, 
we consider the invariants of H of the form: 

^,i(x) > A <^, 2 (x) > 0. 

The following theorem provides a method to determine the invariants of the above form. 

Theorem 2. Let H : (V, L, 7, 0,2? , ^,£q) be a hybrid system. Suppose for each location £ G L, 
there exist two functions (^i(x), <£>£ j2 ( x ) satisfying the following conditions: 

(i) 6 h^o.iW >0A%, 2 (x) >0, 

(ii) <pe,i(x) > 0A^ 2 (x) >0Ag(£,£')Ap(£,£') \= <^' ; i(x') >0A^f i2 (x') > 0, for any transition 

{£,£',9, P) going out of £, 

(iii) <p t ,i(x) = A <pt,i(x) > A *(£) |= <^.i( x ) > 0, 

(iv) <p ltl (x) > A <^. 2 (x) = A *(£) h M x ) > °> 

then </?^i(x) > A ipe^fa) > is an invariant of the hybrid system H at location £. 

Proof. It is easy to prove that Conditions (i), (ii) and (iii) imply that c^i(x) > A ^^(x) > 
satisfies three requirements in Definition 3. □ 
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Similarly, the following corollary shows that the invariant in Theorem 2 can guarantee the 
safety of hybrid systems. 

Corollary 2. Let H : (V, L, T, 6,D, *,£ ) be a hybrid system, and X u {£) = {x G M" : &(x) > 
0} denotes the unsafe region at location £. Suppose for each location £ G L, there exist two 
functions i^i(x), <y3£.2(x) that satisfy the conditions (i-iii) in Theorem 2, and moreover, 

(v) ^,i(x) > A 9^, 2 (x) > h 6(x) < 0, 

i/ien i/ie safety of the system H is guaranteed. 

Remark 1. ^ i(x) > A ^^(x) > in Theorem 2 is an invariant of the hybrid system H at 
location £, but it is not ensured that either tpg i(x) > or ^^(x) > is an invariant of H at the 
location I. This kind of invariants is helpful to look for invariants of the form c\ < /(x) < C2 
where /(x) is a polynomial over R. 

3. Problem Reformulation 

For brevity, we will abuse the notation <pi{x) to represent both the function <p({x) and the 
invariant y«(x) > 0. Clearly, when the functions ^(x) at all locations are identical to y(x), 
then v?( x ) becomes an inductive invariant of the hybrid system. Therefore, in the sequel we only 
discuss how to find invariants v^( x ) for each location fgi, and the problem of computing an 
inductive invariant </?(x) can be handled similarly. Remark that the invariants <pe(x-) or y>(x) are 
also known as barrier certificates in [20]. 

Let us predetermine a template of polynomial invariants with the given degree d. We assume 
that ^(x) = J2 a c aX a , where x Q = x" 1 ■ ■ -x" n , a = (oti, . . . ,a„) G Z" with Y12=i a i d, and 
c a G R are parameters. We can rewrite (pe(x) = c^" • T^(x), where T^(x) is the (column) vector 
of all terms in x\, . . . ,x n with total degree < d, and £ I", with v = is the coefficient 

vector of y^(x). In the sequel, we write y^(x) as y^(x, q) for clarity. 

From Theorem 1, to verify the safety of hybrid system H it suffices to find the invariants (pz(x) 
at each location I G L. The latter problem can be translated into the following problem 

find a e R", W e L 
s.t. e |= <pi (x.,ct ) > o 

<^(x,c<) >0/\g(t,t)Ap(£,H) \=^{x.',c e ,) >0 } (1) 
<^(x, c f )=0A h ^(x, Q) > 
|=^(x,c*) <0 

Suppose that 

6 = {x G R™ : ALi > 0}, X u {£) = {x G R™ : A?=i 0,i( x ) > 0}, 

= {x e R" : ALi iM*) > 0}, ) = {x e R™ : A 4 =i <7«',i(x) > 0}, 

p(£,f)(x,x') = jx'el" : At=iP«',«(x,x') > 0}, 

where £,£' G £, and 6>/(x), <^j(x), ^^(x), <jr^ ,-(x) and pw jlt (x, x') are polynomials over R. 

Let /i(x) G R[x] for 1 < i < m and g(x) G R[x]. Suppose that the set {x G R™ : A™ i /*( x ) > 
0} is compact. According to Stengle's Positivstellensatz, Schmiidgen's Positivstellensatz or Puti- 
nar's Positivstellensatz, if there exist SOS polynomials Oi G R[x] for i = 0, ...,m, such that <?(x) 
can be written as </(x) = cto(x) + 53i=i tJ i( x )/i(x), then the assertion 

m 

A(/ i (x)>0)hff(x)>0 

holds. Therefore, the existence of SOS representations provides a sufficient and necessary con- 
dition of the strict positiveness of <?(x) on the set {x G R" : AI=i/j( x ) — 0}- Moreover, the 
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degree bound of those unknown SOS polynomials <Ji is exponential in n, deg(g) and deg(/j). For 
more details the reader can refer to [13]. Based on the above observation, the problem (1) can be 
transformed into an equivalent SOS programming of the form 

find c e er, Wei 

s.t. tpe {x,c eg ) = cr (x) + J2Li cri(x)6»z(x), 

(fi/(x',ci>) = A«', (x) + J2i=i <W,i(x)s«',i(x) > (2] 

+ Hu=i 7«',«W/0«',„(x,x') + r}M'(x)<Pi(x-,ct), | 
ipt{x.,ci) = 0£,o(x) + J2k=i 0f,fe(x)^,fe(x) + vi(x)tpt(x,Ci) + et,i 
-tpt{x, a) = W,o(x) + Yf J= i Mfj(x)0,j(x) + et,2, 

where ai(x), A#',»(x), 7«',«(x), rjw (x), <^,fc(x), are SOSes in R[x], i^(x) G R[x], and e^i, 

ec,2 £ R+. In practice, to avoid the high computational complexity, we simply set up a truncated 
SOS programming by fixing a priori (much smaller) degree bound 2e, with e £ Z+, of all the 
unknown polynomials. Consequently, the existence of a solution C£ of (2) can guarantee the safety 
property of the given system. 

In the SOS programming (2), the decision variables are the coefficients of all the unknown poly- 
nomials in (2), such as y>f (x, q), <tz(x), j(x). Note that since the coefficients of <pt{x., cg),r]w (x) 
and i^(x) are unknown, some nonlinear terms that are products of these coefficients, will occur 
in the second and third constraints of (2), which yields a non-convex bilinear matrix inequali- 
ties (BMI) problem. In [15], to avoid this BMI problem we strengthened the second and third 
constraints in (1) to 

S (f,<')Ap(<,/)h W (x',c f )>0 and |= ^( x , c t ) > 0, 

respectively, which then results a linear matrix inequality (LMI) problem. For more details, please 
refer to Theorems 5 and 6 in [15]. In this paper, we will discuss in Section 4 how to handle the 
SOS programming (2) directly using the BMI solver or iterative method. 



4. Approximate Solution from BMI Solver 



In Section 3, we have reduced the problem of safety verification of a hybrid system to the SOS 
programming (2) involving BMI constraints. 

Let us first show by an example on how to transform nonlinear parametric polynomial con- 
straints into a BMI problem. 



Example 1. Consider the system x = 2x with location invariant ^ = {x £ 



- 1 < 0}. 



From the discussion in Section 2, to find a polynomial (p(x) satisfying (p(x) > A \& |= (p(x) > 0, 
it suffices to find (p(x) such that 

ip(x) = (fo(x) + 0i(x)(l - x 2 ) + <p 2 (x)ip{x), (3) 

where <j>o(x), <t>i{x), <pz(x) are SOSes. Suppose that deg(</?) = l,dcg(^ ) = 2 and deg(0i) — 
deg(02) = 0, and that if(x) = u + U\X, <f>i = u 2 and <p 2 = v%, with Uo, U\, u 2 , v\ £ R parameters. 
From (3) we have 

4>o{x) — u 2 x 2 + (2 ui - mi vi)x — U2 — Uo vi, 
whose square matrix representation (SMR) [1] is cj>o(x) = Z T QZ 7 where 

i„ 



Q 



—u 2 — Uo V\ U\ 
U\ — \u\ V\ 



2 Ul Vi 

u 2 



and Z 



Since all the <pi(x) are SOSes, we have u 2 > 0, V\ > and Q >z 0, which can be expressed as one 
parametric positive semidefinite matrix 



B(u Q ,u 1 ,u 2 ,Vi] 



u 2 














('1 














-142 - u V\ 


U\ — \u\ V\ 








U\ — \u\ V\ 


u 2 



y o. 
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Therefore, the constraint (3) is translated into a BMI constraint 

2 2 
B = A + Vl Ai + UiA 2+i + ^ UjV! B jA h 0, 

i=0 3=0 

where Ai,Bj^i are constant symmetric matrices. □ 

Similar to Example 1, the SOS programming (2) can be transformed into a BMI problem of 
the form 

inf F(u, v) 

(u,v)eR"+<= ( 4 ) 
s.t. B(u, v) = A + E™i nAi + J2 k J= i VjA m+] + Ei<j< m Ei<j<fe Ui«iB« >= o. 

where Aj, By are constant symmetric matrices, u = . . . , u m ), v = (fi, . . . , fjt) are parame- 
ter coefficients of the SOSes occurring in the original SOS problem, and the objective function 
F(u, v) is a dummy objective function, which is commonly used for optimization problem with 
no objective functions. 

Many methods can be used to solve the BMI problem (4) directly, such as interior-point 
constrained trust region method [14], an augmented Lagrangian strategy [9] and so on. A Matlab 
package PENBMI solver [10], which combines the (exterior) penalty and (interior) barrier method 
with the augmented Lagrangian method, can be applied directly on the BMI program (4). This 
can yield efficiently numerical solutions to the SOS programming (2). 

Alternatively, observing in (4), B(u, v) involves no crossing products like UiUj and ViVj . Taking 
this special form into account, an iterative method can be applied by fixing u and v alternatively, 
which leads to a convex LMI problem [20, 26]. Below is an algorithm. 

Algorithm 1. BMI Solver based on Iteration 

1. [Initialization] Set u = uo. Then (4) becomes to an LMI problem 

inf veRfc F(u ,v) 

s.t. A(u ,v) = A + YhLi UifiAi + Y, k j=i v 3 A rn+J y 0. 

Suppose we obtain a feasible solution v by solving the above LMI problem. 

2. [Fixing v] Find an updated solution u by solving the following LMI problem 

inf ueRm F(u, v) 

s.t. A{u, v) = A + YZLi UiAi + Y?j=i VjA m+J h 0. 

3. [Fixing u] Find an updated solution v by solving the following LMI problem 

inf veRfc F(u, v) 

s.t. A(u, v) = A + Yh=i UiAi + Ylj=l v 3 A m.+j h 0. 

Repeat steps 2 and 3 until a feasible solution (u,v) is found such that B(u, v) >^ 0. 

Remark that, although the convergence of Algorithm 1 can not be guaranteed, the iterative 
method is easier to implement than PENBMI solver. Moreover, from the experiments shown in 
Section 7, Algorithm 1 can yield a feasible solution (u, v) efficiently in practice. 
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5. Exact SOS Recovery 



Xx 









Since the SDP solvers in Matlab are running in fixed precision, the techniques in Section 4 only 
yield numerical solutions of (1). Due to round-off errors, ipe(x, eg) may not be an invariant of 
the given hybrid system at location i, because the constraints in (1) may not hold exactly, as 
illustrated by the following example. 

Example 2. [20, page 31] Consider the following nonlinear system 

x 2 

—X\ + gxf — X2 

we want to verify that all trajectories of the system starting from the initial set will never enter 
the unsafe region X u , where 

9 = {(xx,x 2 ) G R 2 : (xx ~ 1-5) 2 + x 2 < 0.25} 

and 

X u = {(xi,x 2 ) G R 2 : (a* + l) 2 + (x 2 + I) 2 < 0.16}. 

It suffices to find an invariant <p(xi, x 2 ) with rational coefficients, which satisfies all the conditions 
in (1). As stated in [15], we first set up an SDP system using LMI constraints. Apply the SDP 
solver to find a numerical polynomial invariant 

<p(xi,x 2 ) = -1.3686 + 0.62499 x 2 + 1.0669 xxx 2 + 1.5086 x\ - 0.56749 x x x\ 

- 0.15231 x\ - 0.10417 x\ - 0.35564 x\x 2 ~ 0.23739 x\x\ - 0.24152 Xx x\, 

and some associated numerical positive semidefinite matrices . 

Case 1: If we convert the coefficients of ip(x\, x 2 ) to the corresponding rational coefficients sep- 
arately, we obtain 

6843 62499 2 10669 7543 , 56749 9 

ip Xi , X? = 1 Xi H Xi Xo H X? Xi Xo 

Yy ' ' 5000 100000 10000 5000 100000 

15231 , 10417 , 8891 . 23739 2 2 3019 o 

\ x 1 - TTTTTTTTTTT^l _ -— — -Xx X2 - — — — — Xi X2 - ; ^ X^ - 



100000 100000 25000 100000 12500 

However, (^(xx, X2) can not satisfy the conditions in (1) exactly, because there exists a sample 
point p = (— — |) such that the third constraint of (1) can not be satisfied. Therefore, 
(p(xi,X2) is not an exact invariant of this system. 

Case 2: In our former papers [29, 15], we applied Gauss- Newton iteration and rational vector 
recovery to obtain solutions that satisfy the constraints in (1). This technique may fail 
in some cases. Let r = 10~ 2 and the bound of the common denominator of the rational 
polynomial be 100. Then we obtain a rational polynomial 

15 5 2 47 133 2 25 , 

¥>(X1,X 2 ) = -- + -X1 +- Xl X2 + —X2 ~^xix 2 

13 o 9 a 31 o 21 99 21 o 

X2 Xi Xi X2 Xi X9 X1X2 . 

88 88 88 88 88 

However, the subsequent Gauss-Newton iteration and rational vector recovery techniques 
in [29, 15] failed in finding the associated positive semidefinite matrices that satisfy the 
constraints of the polynomial invariant exactly. The reason may lies in that we recover the 
coefficient vector of <^(xx,x 2 ) and the associated positive semidefinite matrices separately. 

In the sequel, we will propose an improved algorithm to compute exact solutions of polyno- 
mial optimization problems with BMI constraints, through a modified Newton refinement and 
rational vector recovery technique applied on the coefficient vector c and the associated positive 
semidefinite matrices simultaneously. 
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5.1. Modified Newton Iteration 



Suppose that applying PENBMI solver or Algorithm 1 in Section 4 yields numerical solutions 
that satisfy (1) approximately. Similar to [7], we now present the modified Newton refinement 
method to refine these solutions. Without loss of generality, we can reduce the problem (1) to 
the following problem 

'find ce W 

s.t. V5i(x,c)>0, (5) 
<£ 3 (x,c) >0 h^2(x,c) >0, 

where the coefficients of the polynomials <^i(x, c), 1 < i < 3 are affine in c. Then, based on SOS 
relaxation, the problem (5) can be further transformed into the following polynomial parametric 
optimization problem 



find cer 

s.t. <Pi(x,c) = mi(x) T • W [1] • mi(x), 

<P2(x, c) = m 2 (x) T • WW . m 2 (x) + (m 3 (x) T • . m 3 (x)) • p 3 (x,c), 
W® >- 0, 1 < % < 3, 



(6) 



involving both LMI and BMI constraints. 

By solving the SDP system (6), we can obtain the numerical vector c and the approximate 
positive semidefinitejiiatrices W™, 1 < i < 3. We first convert to a nearby rational positive 
semidefinite matrix by non-negative truncated PLDL T P T -decomposition, in which all the 
diagonal entries of the corresponding diagonal matrix are preserved to be non-negative. Hereafter, 
denote by 0(x) the rational polynomial m 3 (x) T • Vl/t 3 ! • m 3 (x), and set 9 = ||»"i(x)||| + ||j"2(x)|||, 
the backward error of the numerical solutions of (6), where 



ri(x) = <^i(x,c) - mi(x) T • W [1] • mi(x), 
r 2 (x) = y> 2 (x, c) - m 2 (x) T • ■ m 2 (x) - 

With the numerical c, wWjWl 2 ! and the rational polynomial 
matrix representations in (6) into their SOS forms respectively: 



i(x)<p 3 (x, c). 

(x), we expand the two square 



tpi (x, c) re ^ I ^2p ha x. a J and <p 2 (x, c) 



J2 ( Efc.^ 

P 



(x)^ 3 (x,c), (7) 



where t and k are the ranks of WW and respectively. Apply Gauss-Newton iteration on two 
equations in (7) simultaneously to compute Ac, Api ;tl and Aqj,p such that 



<^(x, c + Ac)w ELl(E Q (Pi,a + Ap,, Q ) x Q ) 2 , 

^ 2 (x,c + Ac)w E; =1 (E fl (?j,/J + A Qj,p)^) 2 + 0(xV 3 (x,c + Ac) 



We update the vector c and matrices by c + Ac and + AW^\ 1 < i < 2, respectively, 
and terminate the Newton iteration when 9 is less than the given tolerance r. In doing so, we 
will obtain the refined solution, the vector c and matrices W^- 2 \ of (6) such that 



pi(x,c) - mi(x) T • . mi(x) re 0, 

<^ 2 (x,c) - m 2 (x) T • ■ m 2 (x) - ^ 3 (x,c)^(x) 

WW 



w 



>Z 0, W T = W, 



w&_ 

and the backward error of the numerical solution satisfies 9 < r. 



(8) 
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5.2. Exact Recovery 

In this section, we will present two error-free algorithms to obtain a rational vector c, which 
satisfies the constraints in (5) exactly. 

Discussed in Section 5.1, we can obtain a refined solution c 6 R fc by applying the modified 
Gauss-Newton iteration. Instead of recovering the entries of c separately, we can deploy simul- 
taneous recovery technique [12] to achieve a rational vector c near to c for a given bound D of 
the common denominator of c. Then we need verify whether c is an exact solution of (5), that 
is, whether c satisfies 

<£i(x, c) > and ^ 3 (x, c) > |= (£> 2 (x, c) > 0, 

or equivalently, both 

<f i (x, c) < and (^3 (x, c) > A <^ 2 (x, c) < 

has no real solutions. Clearly, verifying the solution c is equivalent to determining two con- 
stant semi-algebraic systems have no real solutions, which can be verified by Maple packages 
Regular Chains, DISCOVERER [30] and RAGLib [3]. 

In this paper, we focus on retrieving exact SOS representations for (6). We will discuss how to 
recover from c, W M , , the rational vector c and rational positive semidefinite matrices 
and that satisfy exactly 

< y si(x,c)-mi(x) T - W [1] -mi(x) = and p 2 (x, c)-m 2 (x) T - W [2] ■m 2 (x)-^ 3 (x, c)</.(x) = 0. (9) 

Since the equations in (9) are affine in entries of c and W\,W 2 , one can define an affine linear 
hyperplanc 

C = |c, W [1] , W [2] I 951 (x,c) - mi(x) T • W [1] -mi(x) =0, 

^ 2 (x, c) - m 2 (x) T • • m 2 (x) - p 3 (x, c)0(x) = } . (10) 

Note that the hyperplane (10) can be constructed from a linear system Ay = b, where y consists 
of the entries of c, W^. If A has full row rank, such a hyperplane is guaranteed to exist. 
Then the rationalized SOS solutions of (6) can be computed by orthogonal projection if the matrix 
W in (8) has full rank, or by the rational vector recovery otherwise. 
Case 1: W is of full rank 

Suppose the refined matrix W in (8) is of full rank numerically, namely, the minimal eigenvalue 
of W is greater than the given tolerance r. In this case, for a given bound D of the common 
denominator, we apply orthogonal projection technique to obtain the rational vector c and rational 
matrix W, which lie on the affine linear hyperplane defined by (10). This projection can be 
achieved by exactly solving the least squares problem: 

nun ||c-c||!+ IIW-WII! 

c, W 

s. t. ¥>i(x,c) = mi (x) T • WW ■mi(x), f ( U ) 

y 2 (x,c) = m 2 (x) T • W [2] ■ m 2 (x) + </>(x)^ 3 (x, c). 

For the rational solution c and W of (11), we compute the exact PLDL T P T -decomposition to 
check whether W is positive semidefinite. If so, then c is verified solution of (6). 

Theorem 3. Let c, W be the refined numerical solution of ( 6) with the backward error 9 < r . 
and Ay = b be the linear system associated to the hyperplane (10). Suppose that c and W are the 
optimal rational solutions of the least square problem (11). Let A £ K>o be the minimal eigenvalue 
of W . If A has full row rank and A > 2?7k 2 (A)t 2 with 77 = |]c||| + ||W|||., then W is positive 
semidefinite, and c is a certified solution of (6). 



10 



Proof. Clearly, c and W satisfy the equations in (9). Let z and z be the vectors consisting 
of the entries in c, W and those in c, W, respectively. Since A has full row rank, we have 
\\Az — &H2 = < r, and then Az = b. According to the perturbation result in [4, Th.5.7.1] for 
full rank underdetermined systems, we have 

||z-z|| <MA)r)||z|| 2 +0(r 2 ). (12) 

From (12) and the assumption A > 2 ^k 2 {A) r 2 , we have \\W — W\\ F < 2 ||z — z||| < A, where the 
last inequality follows from that 0(t 2 ) is negligible in comparison with r\ when r is very small. 
Let A be the minimal eigenvalue of W. By Wielandt-Hoffman theorem [4, page 395], we have 
|A - A| < \\W - W\\% < A, which concludes that W h 0. □ 

Remark 2. For a numerical solution of (6), 2t7k 2 (A)t 2 in Theorem 3 is a constant. Therefore, if A 
has full row rank and the given tolerance r is small enough, then the orthogonal projection method 
by solving exactly least squares problem (11) will always find a rational positive semidefinite 
matrix W, therefore yielding a certified rational solution c of (6). 

Case 2: W is singular 

It is pointed out in [8] that, the orthogonal projection by solving (11) may fail if the resulted 
rational matrix W satisfying (9) is positive semidefinite but not strictly positive definite. In this 
case, we need to explore rational vector recovering technique to obtain the rational c and W. 
Unlike in [15], the technique to be presented will recover the rational vector c and matrix W 
simultaneously. More specifically, given a common denominator bound of the entries in c and 
W, we employ the simultaneous Diophantinc approximation algorithm [12] for c and the singular 
matrices in {W' 1 ', W' 2 '}. There are two cases to be addressed. 

Case 2.1: Both and are singular. We apply rational vector recovery directly on c and 
W simultaneously to obtain c and W. 

Case 2.2: Only one matrix, say, is singular. We apply vector recovery technique on c and 
simultaneously, and then check whether >z by PLDL T P T -decomposition. If so, 
the rational matrix can be obtained by orthogonal projection by solving the following 
exactly least squares problem: 

min \\WW -W^\\j, ) 

wm > 

s.t. <^ 2 (x,c) = m 2 (x) T • W [2] ■ m 2 (x) + </? 3 (x, c)0(x). J 

Remark 3. Similar to (6), we first apply the modified Newton iteration to compute the refined 
numerical solutions c and WW of (2) respectively. For the given tolerance and common denom- 
inator, we employ rational vector recovery technique simultaneously for c and all the singular 
matrices among the to find rational c and the associated rational positive semidefinite ma- 
trices. Finally, orthogonal projection is applied on the remaining nonsingular matrices to obtain 
the rational non-singular positive definite matrices. 

5.3. Algorithm 

The results in Sections 5.1 and 5.2 yield an algorithm to find exact solutions to (6). 

Algorithm 2. Verified Parametric Optimization Solver 

Input: ^ a polynomial optimization problem of the form (6). 

^ fle Z>o: the bound of the common denominator. 

^ e e Z>o: the degree bound 2e of the SOSes used to construct the SOS program- 
ming. 

^ t£ K>o : the given tolerance. 
Output: ^ the verified solution c of (6) with the 1 < i < 3 positive semidefinite. 
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1. [Compute the numerical solutions] Set up the SOS programming (6) with the degree 
bound 2e and apply PENBMI solver or Algorithm 1 to compute its numerical solutions. If 
the system has no feasible solutions, return "we can't find solutions of (6) with the given 
degree bound 2e". Otherwise, obtain c and W^' fc 0, 1 < i < 3. 

2. [Compute the verified solution c] 

(2.1) Convert to a nearby rational positive semidcfinitc matrix W' 3 ' by non-negative 
truncated PLDL T P T -dccomposition. 

(2.2) For the tolerance r, apply the modified Newton iteration to refine c and W^. 

(2.3) Determine the singularity of and with respect to r. 

Case 1: Both and are of full rank. The rational vector c, and rational matrices 
and can be obtained by orthogonal projection for a given bound D of 

the common denominator Check if >r for i = 1,2. If so, return c, and 
1 < i < 3. Otherwise, return "we can't find the solution of (6) with the given 

degree bound" . 

Case 2: Both W [1] and are singular. Obtain c and as described in Case 2.1 of 
Section 5.2. 

Case 3: Only one of and WW is singular. Obtain c and as described in Case 2.2 
of Section 5.2. 

6. Computational Issue 

In practice, it may happen that the unsafe region is too big to find an appropriate invariant 
for verifying the safety of the given system. To deal with this issue, we can divide the original 
unsafe region into two parts, which can be achieved by bisection method through computing the 
minimum and maximum values of some variable xi. Then we compute the invariants that verify 
safety for two unsafe sub-regions X Ut i,X U)2 , respectively. This procedure can be easily repeated 
until the safety property have been verified for all of those unsafe sub-regions. The process for 
splitting the initial set O can be handled similarly. 

The following example is presented to illustrate the above technique. 

Example 3. [22, Example 2] Consider the following nonlinear system 



x\ 




Xi 


- x 2 


X 2 




Xi 


+ x 2 



within the region * = {(xi,x 2 ) S K 2 : < x 1 < 4 A < x 2 < 4}. We want to verify that all 
trajectories of the system starting from the initial set O will never enter the unsafe region X u , 
where 6 = {{x x ,x 2 ) £ M 2 : 2.5 < x 1 < 3 A x 2 = 0} and X u = {(x u x 2 ) e R 2 : x\ < 2}. 

For X u , PENBMI solver or Algorithm 1 can not yield a feasible solution of the associated SOS 
program. Therefore, it is impossible to obtain the invariant to verify the safety property for this 
system. Here we split the original unsafe region into X U: i = {(xi,x 2 ) 6 K 2 : x\ < 2 A x 2 < 0} 
and X Uy2 = {(xi,x 2 ) 6 R 2 : x\ < 2 A — x 2 < 0}. Then, our certified method, combining modified 
Gauss-Newton refinement with exact recovery technique, can yield the verified invariants with 
degree 2 for X Ui i and X u ^ 2 , respectively. Thus the safety of this system is guaranteed. 

7. Experiments 

Let us present some examples of safety verification for hybrid systems. 

Example 4 (Example 2 Revisited). Consider the safety verification problem of Example 2. 
We now apply Algorithm 2 to compute two invariants (p\ and <p 2 , that are subject to the strength- 
ened LMI constraints in [15] and the BMI constraints in (2), respectively. 
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According to Theorem 6 in [15], we construct an SOS system involving only LMI constraints. 
Applying the improved Gauss-Newton iteration and rational vector recovery techniques, we find 
an exact invariant (/?i(xi,X2) > until deg(<^>i) = 4, and the corresponding SOSes. Here we only 
list the invariant 

53 8 , 59 2 2 , 4 , 

(jSl Xl , Xo ) — 1 XT H XT Xn XT 

V i> ^ 39 13 1 39 2 13 2 39 1 

14 22 2 14 3 3 2 2 3 , 

+ YE Xl X2 ~39 Xl X2 ~ 39 Xl X2 ~ n Xl X2 ~ 13 Xl X2 ' 

The invariant </3i(xi, x%) > guarantees the safety of the given system. 

Alternatively, from Theorem 1 we construct another SOS system with BMI constraints and 
then apply Algorithm 2 to obtain an exact invariant with degree 2: 

151 62 152 106 4 2 

(fio ( Xl , Xo ) = Xo Xl Xl Xo Xi . 

' ' 99 33 99 99 9 1 

This proves the safety of the given system. 

Consider again the system in Example 2 with the same initial set O but a larger unsafe region 
X u = {(xi,x 2 ) € R 2 : (xi + l) 2 + (x 2 + l) 2 < 1}. Similarly, we construct the SOS systems with 
LMI constraint and BMI constraint respectively, and apply Algorithm 2 to compute the invariants 
£>i(xi,.X2) > with deg((J5i) = 6, and ^2(^1,^2) > with deg(£>2) = 4, where 



1714 1856 3 6160 3 3927 2 1507 3 75 5 426 

3209 + 3209 Xa + 3209 Xl ~ 3209 X2 ~ 3209~ Xl X2 + 3209 X2 + 3209 



16 73 71 147 , 299 2 2 347 , 211 4 

g9 2 xi,x 2 ) = Xi H xo x, + ■ • • 1, 1, xi x 2 x, . □ 

YK ' ; 25 50 50 100 1 100 1 2 100 2 100 1 

From Example 4, we see that the invariants obtained from the BMI constraints are of lower 
degree than those obtained from the LMI constraints. 

Example 5. Consider a hybrid system [20] depicted in Figure 1, where 



/i(x,d) = 



•' : 2 

-Xi + X 3 

xi + (2x 2 + 3x 3 )(l + x 2 ) + d 



/a(x,d) = 



X2 

-Xl + x 3 
-Xi — 2x2 — 3x 3 



with parameter d satisfies — 1 < d < 1. The system starts in location t\, with an initial state in 

9 = {(xi,x 2 ,x 3 ) e M 3 : x 2 + Xj + x 2 < 0.01}. 

Our task is to verify the system never reach the states of 

X u {i 2 ) = {(xi,x 2 ,x 3 ) G R 3 : 5 < xi < 5.1 V -5.1 < Xj < -5}. 

To prove the safety of the hybrid system, it suffices to find the corresponding invariant polyno- 
mials v ? i( x ) and ^(x) at locations ii and £2, which satisfy all the conditions in Theorem 1. 
Applying Algorithm 2, we obtain the invariant polynomials with rational coefficients 

_ . , 53 39 5 1 1 2 3 2 3 2 

05l(x)= XlX 2 H X1X3 X2X3 Xl X2 X3 , 

y y ' 44 88 88 88 44 88 88 ' 

„ , . 129 1 1 1 1 2 1 2 1 2 

W9<x = X1X2 H X1X3 -I X2X3 -I xi H X2 -I X3 . 

y - y ' 22 88 88 88 88 88 88 

Moreover, £>i(x) and ^i(x) satisfy the conditions in Theorem 1 exactly. Therefore, the invariants 
can guarantee the safety of the hybrid system. 
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0.99 <x\ + O.Ola;! + 0.01.x| < 1.01 




0.03 < x\ + x\ + x\ < 0.05 



Figure 1 : Hybrid system of Example 3 

Table 1 shows the performance of Algorithm 2, for safety verification of some interesting 
benchmark examples. All the computations have been performed on an Intel Core 2 Duo 2.0 
GHz processor with 2GB of memory. Examples 1-3 correspond to Examples CLOCK, FOCUS, 
ECO in [23], and Examples 4-7 correspond to [21, Example 3], [32, Example 3], [22, Example 
11] and [31, Example 2]. All these examples except Example 2 are nonlinear systems. In all the 
associated SOS programmings, the degree bound of SOSes is e = 6 and we set D = 1000 and 
r = 10~ 10 in Algorithm 2. Here \£\ and n denote the number of locations and the number of system 
variables respectively; BMI Solver refers to the method used to obtain the numerical solutions of 
the given BMI problems; Num. Para, is the number of decision variables appearing in the BMI 
problem; deg((^) and deg(^LMi) denote the degrees of invariants computed by Algorithm 2 and 
by the algorithm used in [15] with LMI constraints, respectively. Fail means that the algorithm 
in [15] fails to find invariants with degree < 6. Time is that for the entire computation run 
Algorithm 2 in seconds. 



Ex. 


14 


n 


BMI solver 


Num. Para. 


dcg(v?) 


deg(v3 L Mi) 


Time (s) 


1 


1 


2 


PENBMI 


30 


2 


3 


6.17 


2 


1 


2 


Alg. 1 


16 


2 


Fail 


4.56 


3 


2 


2 


Alg. 1 


70 


2 


2 


14.21 


4 


1 


2 


PENBMI 


18 


2 


Fail 


4.31 


5 


1 


2 


PENBMI 


41 


3 


Fail 


7.99 


6 


1 


3 


Alg. 1 


21 


2 


2 


5.45 


7 


1 


2 


Alg. 1 


24 


2 


2 


5.62 



Table 1: Algorithm Performance on Benchmarks 



8. Conclusion 

In this paper, we present a symbolic-numeric method on safety verification of hybrid systems. A 
numerical invariant of a hybrid system can be obtained by solving a bilinear SOS programming 
via PENBMI solver or iterative method. Then a method based on modified Newton iteration and 
rational vector recovery techniques is deployed to obtain exact polynomial invariantswith rational 
number coefficients. Some experimental results are given to show the efficiency of our method. 
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